VIBRATIONAL RELAXATION OF DIATOMIC MOLECULES 
IN SOLIDS AT LOW TEMPERATURES 


LAWRENCE L. HALCOMB 
AND 

DENNIS J. DIESTLER 

DEPARTMENT OF CHEMISTRY 
PURDUE UNIVERSITY 

WEST LAFAYETTE, INDIANA 



VIBRATIONAL RELAXATION OF DIATOMIC MOLECULES 


IN SOLIDS AT LOW TEMPERATURES 

* 

Lawrence L. Halcomb and Dennis J. Diestler 
Department of Chemistry 
Purdue University 
West Lafayette, Indiana 47907 


Prepared for delivery at the joint NASA/Goddard-CDC Symposium 
on CYBER 205 Applications, Lanham, Maryland, October 11-12, 1983. 


Abstract 

A miscroscopic dynamical treatment of chemical systems comprising both 
light particles that require a quantal description and heavy ones that may be 
described adequately by classical mechanics has recently been presented 
[J. Chem. Phys. 78, 2240 (1983)]. The application of this ' ' hemiquanta 1 ' ' 

method to the specific problem of the vibrational relaxation of a diatomic 
molecule embedded in a one— dimens ional lattice is presented. The vectorization 
of a CYBER 205 algorithm which integrates the 10 3 -10 4 simultaneous 
' ' hemiquanta 1 ' ' differential equations is examined with comments on opti- 
mization. Results of the simulations are briefly discussed. 
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I. Introduction 


A microscopic dynamical description of a chemical system composed of both 
light particles that require a quantal description and heavy ones that may be 
described adequately by classical mechanics has been proposed recently [J. 
Chem. Phys. 78 , 2240 (1983)]. The description consists of a self-consistent 
set of ' * hemiquantal ' * equations (HQE) arrived at by taking a partial classical 
limit of Heisenberg's equations of motion for the system. In form, the HQE 
appear to consist of Heisenberg's equations for the light particles coupled to 
Hamilton's equations for the heavy particles. The coupling is self-consistent 
in that there is an instantaneous feedback between the light and heavy 
subsystems, with total energy and probability of presence of the quantal 
subsystem being conserved. 

This paper will focus on the numerical solution of the HQE on the CYBER 205 
for the special case of a diatomic molecule embedded in a cold, one-dimensional 
lattice. In Section II, we detail the model and specific form of the HQE, 
while the CYBER 205 algorithm and steps taken to optimize performance are 
included in Section III. Results of the simulations and some discussion of 
their physical significance are presented in Section IV. 

II. Model and Equations of Motion 

Figure 1 depicts the physical situation, i.e. a single diatomic molecule BC 
occupying a substitutional site in an otherwise pure one-dimensional lattice of 
atoms Aj the end atoms of the lattice are assumed free. So that the normal 
modes of the lattice are known analytically, the mass of BC is taken to be 
equal to that of A. The heavy, classically behaving degrees of freedom are 
considered to be the displacements (u^) of the lattice atoms, including the 
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center of mass of BC, from their equilibrium positions. The internal vibration 
(q) of BC is treated qnantally and, for simplicity, as a harmonic, two-state 
system. Ve assume that only nearest-neighbor atoms interact with one another: 
A-A interactions are harmonic ^ A-B and A-^C interactions are approximated by 
Morse potentials. 

Under these conditions, the HQE take the form 

c.(t) = -ilfhe.c.U) + 2 V ij <(tt k (t)})o j (t)] 

j 

u (t) = p (t)/m (1) 

l l A 

Pf (t) = ” f„ 0({u.(t))) + ) o.»(t)o. (t)F. .. (tu <t)}) . 
i ou. j L j x ljk m 

jk 


Here c^ is the occupation probability amplitude for quantal state i^ p^ is 
momentum conjugate to a.J U is the harmonic part of the potential, i.e. 


the 


n-2 


N-l 


°"l l l ( " i+ i - V 2 + } 

i=l i*n+l 

where N is the number of lattice atoms. F is the quantal force 
defined by 


(2) 


F. .. = 3V. . / 3u. 

ijk ij k 


( 3 ) 


where 
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(4) 


V<V> “ <ilV AB + V AC U> ’ 


and the Morse potential V An is explicitly 

Ad 


V AB = D AB teIp[ ‘ a AB ( Wl + L ‘ Y B q,1 ' 1) 


(5) 


with a similar expression for V Ar . 

Since the c^ are complex, the HQE consist of 2N+4 coupled first-order 
ordinary differential equations. Given initial conditions appropriate to the 
physical situation, we can integrate these numerically by standard techniques. 
Our principal problem now is to develop and optimize an algorithm appropriate 
to the CYBER 205. 


III. CYBER 205 Algorithm 


The HQE [Eqs. (1)] can be cast in terms of the vector differential equation 
X = f(X(t)) , defined by 

x (t) = f (x , ...» x ), x (0) = x° , 

1 11 n l 1 


x (t) = f (x , . 
n n l 


. . , x ) , x (0) = x 
n n n 


The vector X can be written as 


(6) 


X = [C,U,P] where, for example, 

C = [Cl, C2 , C3 , C4J . (7) 
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From experience, we have found the HOE extremely well-behaved. Therefore, they 
can be handled with a relatively simple differential equation solver. We 
employ the familiar fourth-order Runge-Kutta algorithm (RK4) which, for our 
case, is summarized by the following equations: 

= T f(X) 

K 2 = T f(X + K 1 /2) 

X 3 = T f(X + K 2 /2) (8) 

X 4 = T f(X + X 3 ) 

X[ (n+1 ) T] = X(nT) + (X^K^/6 + (* 2 +K 3 )/3 

where T is an appropriately chosen time step. Our choice of RK4 is guided by 
several cons iderations \ it is quite stable, self-starting and easily coded for 
the CYBER 205. In addition, we need no direct method of estimating truncation 
error since we can calculate total energy and probability of the system as a 
check. Eventually, the RK4 algorithm will be used to calculate input values 
for a more sophisticated predictor-corrector routine. 

Since our simulations require widely varying amounts of memory, we would 
like to assign storage at execution time. Clearly, the vector pipelines are 
used more efficiently if the entire derivative vector is manipulated at once. 
If we are to deal almost exclusively on the dynamic stack, we need a method of 
parsing the vector X into subvectors C,U,P which can then be handled 
independently. This ''breaking up'' is accomplished by building descriptors 
using SHIFT and OR operations on an integer equivalenced to a descriptor which 
points to an area in dynamic space. The subroutine BREAKUP is presented in the 
Appendix. This routine allows the RK4 mainline to allocate storage dynamically 
while permitting the derivative routine to access each subvector individually. 
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We now concentrate on the vector function subprogram that calculates the 
derivative f(X). In our model, the four probability amplitudes must be 
accessed individually each time the function is called. Rather than waste a 
vector instruction to store the subvector C in a temporary array, it is faster 
and more convenient to use the following sequence of hardware calls to load 
them directly into registers: 

ASSIGN TEMP, C 

CALL Q8L0D (TEMP, , Cl) 

CALL Q8IX(TEMP, 64) 

CALL Q8L0D( TEMP,, C2) , etc. 

The constants needed to calculate the potential and force functions are 

computed in advance and passed via labeled common. By reviewing an assembly 

listing of the program, one can minimize the number of loads necessary to 

access these constants. The evaluation of U is easily done by a vector 

multiplication with a stored reciprocal mass. 

♦ 

P can be conveniently calculated by evaluating the derivative of a fully 
harmonic potential U' . Thus we have 

- U'({u.}) = k(-2u.+u. ,+u.. .) where 

ou. 1 1 l-l 1+1 

i J 

u o ■ U 1 ■ Vi " V <9) 

which can be effected by two vector additions and two vector multiplications as 
follows : 
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~To u ' (fu j }) = dtemp(i;n) - k«(-2.*utemp(i;n> + utemp(o;n) + utemp(2;n>) 

» 

where UTEMP is a temporary array set to the current values of U. Finally, P is 
obtained by replacing the n-1, n, and n+1 elements of UTEMP by the proper 
values reflecting the Morse potentials at the diatomic. To accomplish this, it 
is necessary to access the five displacements {u^, i = n-2, n+2} . Alternative- 
ly, descriptors could be built to define the necessary vectors on U and the 
values stored in UTEMP. In this case, hardware calls would be required to set 
the first and last elements of UTEMP, to access the five elements of U around 

u , and to store values in the three middle positions, 
n 

The conservation of total energy and probability gives us two necessary 
criteria to check the accuracy of the numerical solution. The total energy is 
given by 


E = U( {u. }) + P* P/ ( 2m ) 
1 A 


+ I c 1 2 e + j c I 2 e 
o o 11 


+ ,c J 2v ™ + 2Re{c •c 1 )V lft + 

o oo o 1 10 1 11 


( 10 ) 


while total probability is simply 

P = ic I 2 + I c 1 1 2 , (11) 

o 1 

which must remain unity. These checks were made every 1000 iterations using 
values calculated in the first pass through the derivative routine. To 
calculate U' [Eq. (9)], the following code is used: 
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ASSIGN TEMP# .DIN. N-l 


TEMP= Q8VDELT(U;TEMP) 

EU= (K/2) * Q8SD0T( TEMP, TEMP) . 

In Table I, sample iteration times and estimates of floating point 
operations per second are given. The timings are for loops without I/O or 
accuracy checks. The results of several simulations are presented in the next 
Section. 

IV. Results of Simulations 

Our simulations all take the diatomic to be in its excited state and the 

lattice to be at OK initially. This means that all elements of X(0) are zero# 

except the real component of 0^(0), which is unity. The time step size is .01 

a) ^ , where m is the transition frequency of the diatomic. The quantity of 

2 

principal interest here is |c^| » the probability of the diatomic being 
excited. The physical constants for the system# which are chosen roughly to 
mimic HC1 in Ar, are listed in Table 2. The only variable quantities are ti> and 
N. The transition frequency is chosen low in order to observe relaxation on 
the time-scale of the simulation. 

2 

Figure 1 displays plots of Jc^f versus time for a sampling of simulations. 
Frames (a)-(c) demonstrate the effect of increasing the diatomic's transition 
frequency w, (cm *) holding the number of lattice atoms fixed. It appears that 
the rate of loss of energy from the diatomic increases with increasing 
frequency up to a point. In fact, frame (c) suggests that the diatomic evolves 
to a metastable state in which it loses no further energy. To test this 
hypothesis, we increased the number of lattice atoms to N = 2000. The result, 
shown in frame (f), bears this notion out. For purposes of comparison, we 
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include a simulation for a smaller lattice (N = 200) . Here we see the effect 
of a pulse, which bounces back and forth, interfering with the monotonic 
relaxation of the diatomic. 

V. Conclusion 

These simulations represent the first application of a new description of 

the dynamics of chemical processes. Most previous approaches employ long-time 

asymptotic approximations, in which the coupling between the subsystems is weak 

and the decay is therefore very slow on the time scale of molecular motions 

(10 14 s) . The advancement of ultrafast laser spectroscopy now allows chemists 

—12 

to monitor directly fast relaxation processes (10 s) . In this regime, the 
coupling is more significant, and accurately solving the equations of motion 
becomes crucial. The HQE can be used for this purpose. However, any practi-cal 
implementation will require a vector processor, such as the CYBER 205. Our 
calculations would be essentially impossible on Purdue University's 

6500/6500/6600 system, for example. The calculations would take 50-100 times 
longer, even if the storage for the vectors were available. 

The main feature of our CYBER 205 algorithm is a mainline that assigns 
storage at execution time. The vector function subprogram that evaluates the 
derivative can access the subvectors individually while the mainline processes 
the entire vector. This is accomplished by building the appropriate 
descriptors using the BREAKUP subroutine (see Appendix). 

Some preliminary results were presented in Section IV. Future research 
will deal with the actual mechanism of energy exchange between the two 
subsystems. Also planned are some N-state models with applications in surface 
chemistry . 
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Table 1. Increase of calculation speed with increase 
of number of equations 




Equations 

Iteration Time 

Mega FLOPS 


24 

.157 

ms 

6.1 

204 

.204 

ms 

22.8 

804 

.256 

ms 

37.9 

2004 

.671 

ms 

69.3 

4004 

1.19 

ms 

77.7 

10004 

2.75 

ms 

83.8 

20003 

5.37 

ms 

85 .6 
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Table II. Parameters of model system 


= 9.25 x 10 ^ ergs D^ c 

a AB = 1 - 83 1 10 * cm_1 a AC 

2 

k = 814 ergs /cm m^ 

m B = 1.67 x 10" 24 g m c 
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Appendix 

SUBROUTINE BREAKUP (X, NSUB , LENSUB, DESSUB , NDIM) 

IMPLICIT INTEGER ( A- Z) 

C 

C BREAKUP- TAKES A DESCRIPTOR (X) AND MANUFACTURES OTHER 
C DESCRIPTORS [DESSUB (N)] THAT POINT TO SUBVECTORS OF 

C LENGTHS LENSUB (N) WHICH COMPRISE THE VECTOR POINTED 

C TO BY X 

C 

C ARGUMENTS : 

C X- DESCRIPTOR TO BE 'BROKEN UP' 

C NSUB- NUMBER OF SUBVECTORS 

C LENSUB- ARRAY CONTAINING THE SUBVECTOR LENGTHS 

C DESSUB- ARRAY CONTAINING THE RESULTING DESCRIPTORS 

C NDIM- DIMENSION OF LENSUB AND DESSUB 

C 

DESCRIPTOR D,X,DESSUB(NDIM) 

DIMENSION LENSUB (NDIM) 

EQUIVALENCE ( D , DTEMP ) 

ASSIGN D,X 

ADD= SHIFT ( SHIFT ( DTEMP, 16 ), -16) 

DO 100 N=1,NSUB 

LENGTH^ SHIFT ( LENSUB (N), 4 8 ) 

DTEMP= 0R( ADD, LENGTH ) 

ASSIGN DESSUB (N),D 
A DD= ADD + 64*LENSUB(N) 

100 CONTINUE 

RETURN 

END 
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